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Abstract For many systems of differential equations modeling problems 
in science and engineering, there are natural splittings of the right hand 
side into two parts, one non-stiff or mildly stiff, and the other one stiff. 
For such systems implicit-explicit (IMEX) integration combines an explicit 
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Q scheme for the non-stiff part with an implicit scheme for the stiff part. In a 

recent series of papers two of the authors (Sandu and Zhang) have developed 

IMEX GLMs, a family of implicit-explicit schemes based on general linear 

methods. It has been shown that, due to their high stage order, IMEX 

C — GLMs require no additional coupling order conditions, and are not marred 

^ by order reduction. This work develops a new extrapolation-based approach 

• to construct practical IMEX GLM pairs of high order. We look for methods 

(^ with large absolute stability region, assuming that the implicit part of the 

^ method is A- or L-stable. We provide examples of IMEX GLMs with optimal 

stability properties. Their application to a two dimensional test problem 

• ^ confirms the theoretical findings. 
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1 Introduction 

Many practical problems in science and engineering are modeled by large sys- 
tems of ordinary differential equations (ODEs) which arise from discretization 
in space of partial differential equations (PDEs) by finite difference methods, 
finite elements or finite volume methods, or pseudospectral methods. For 
such systems there are often natural splittings of the right hand sides of the 
differential systems into two parts, one of which is non-stiff or mildly stiff, 
and suitable for explicit time integration, and the other part is stiff, and 
suitable for implicit time integration. Such systems can be written in the 
form 

' y'it) = f{y{t))+g{y{t)), t e [to,T], ^^ ^^ 

y{to) = z/o, 

where f{y) represents the non-stiff processes, for example advection, and 
g{y) represents stiff processes, for example diffusion or chemical reaction, in 
semi-discretization of advection-diffusion- reaction equations [18] . 

Implicit-explicit (IMEX) integration approach discretizes the non-stiff 
part f{y) is with an explicit method, and the stiff part g{y) with an implicit, 
stable method. This strategy seeks to ensure the numerical stability of the 



solution of (1.1 ) while reducing the amount of implicitness, and therefore the 
overall computational effort. IMEX multistep methods were introduced by 
Crouzeix [13j and Varah |23] and further analyzed in [21 [E]. IMEX Runge- 
Kutta methods have been investigated in [H [TOl 1201 [221 ESI EO] • 

In a recent series of papers the last two authors and their collaborators 
have proposed the new IMEX GLM family of implicit-explicit schemes based 
on general linear methods. A general formalism for partitioned GLMs and 
their order conditions was developed by Zhang and Sandu [28j. The par- 
titioned method formalism was then used to construct IMEX GLMs. The 
starting and ending procedures, linear stability, and stiff convergence prop- 
erties of the new family have been analyzed. Zhang and Sandu examined 
practical methods of second order in [2^ and of third order in [28] . A class 
of IMEX two step Runge-Kutta (TSRK) methods was proposed by Zharovski 
and Sandu [22] • 

The results in f27\ [251 [2S] prove that the general linear framework is 
well suited for the construction of multi-methods. Specifically, owing to the 
high stage orders, no coupling conditions are needed to ensure the order of 
accuracy of the partitioned GLM [21] . In addition, it has been shown that 



IMEX GLMs are particularly attractive for solving stiff problems, where 
other multistage methods may suffer from order reduction [28] . 

This paper extends our previous work [271 128| 129] and develops a new 
extrapolation-based approach for the construction of practical IMEX GLM 
schemes of high order and high stage order. 

The organization of this paper is as follows. General linear methods and 
the implicit-explicit variants are reviewed in Section |2} The new extrapolation- 
based IMEX GLMs are derived in Section [3} and their order conditions are 
presented. The stability analysis is performed in Section |4] and specific meth- 
ods are constructed in Section |5} Numerical experiments are presented in 
Section |6| and Section [7] gives some concluding remarks and plans for future 
work. 



2 Implicit-explicit general linear methods 

In this section we briefly review GLMs and the IMEX GLM family. 

The GLMs for ODEs were introduced by Burrage and Butcher [1] and 
further investigated in [3l [5l [H [HI [T2|, [151 HSl [E] ■ We also refer the reader to 
the review article [8] and the recent monograph [19] and references therein. 



A diagonally implicit GLM for (1.1) is defined by 



i=i i=i 

i=i j=i 

(2.1) 

n = 0, 1, . . . , A^— 1. Here, iV is a positive integer, h = (T—to)/N, tn = tQ+nh, 

n = 0, 1, . . . , A^, Fj are approximations of stage order q to y(t„ + q/i), 

i.e., 

yj"^'^ = y{tn + Q/i) + 0(/i«+i), 2 = 1, 2, . . . , s, (2.2) 

yl' are approximations of order p to the linear combinations of the derivatives 
of the solution y at the point t„, i.e., 

yf^ = Y.q^kh^y^'\tn) + 0{h^^^), ^ = l,2,...,r, (2.3) 

k=0 



and y is the solution to (1.1). These methods can be characterized by the 



abscissa vector c = [ci, . . . , c^] , the coefficient matrices A = [aij] 



U 



[aij\ e 



B 



V 



[aij\ e 



the vectors 



qo, • . . , Qs € M*" defined by q^ = [(lj/i\i<j<r^ and four integers: the order p, the 
stage order g, the number of external approximations r, and the number of 
stages or internal approximations s. 



The method (2.1) can be written in a compact form 






,["+!] 



(2.4) 



n = 0,l, 



A^ — 1, and the relation (2.3) takes the form 



(2.5) 



Applying (2.1) to the basic test equation y'{t) = Xy{t), t > 0, A G C, 
leads to the recurrence equation 

y["+i] = S(%M^ n = 0,l,..., 
z = hX, with the stability matrix given by 

S{z) =Y + zB{I-zA)-^V. 
We also define the stability polynomial ri{w, z) by 

Tjiw, z) = det (wl — S{z)) . 



(2.6) 



(2.7) 



The region of absolute stability of the method (2.1) is the subset of the 
complex plane 

A = {z ^ C : all roots Wi{z) of ri{w, z) are in the unit circle}. (2.8) 

The traditional concepts of A(a)-stability, A-stability, and L-stability apply 



directly to GLMs via (2.8). 



In this paper we will examine only methods of high stage order, i.e., 
methods where q = p — 1 or q = p. It has been shown in [HI El [H] that the 



GLM (2.1 ) has order p and stage order q = poTq = p — lii and only if 



k\ 



J2ji^k-i-kBc'-'-k\Vcik 



£=0 



0. 



0. 



fc = 0, l,...,g, and (2.9) 



fc = 0,l,...,p. 



(2.10) 



An IMEX-GLM [291 Definition 4] has the form 



y[n+i] ^ /i(A<=^P I) /(F["+il) + /i(Ai"^P 1) ^(F["+il) + (U ® I)y["] 



/i(B<=^P ® I) /(F["+il) + /i(B''^P ® I) c/(r ["+!]) + (V ® I)|/W^ 

(2.11) 
B™P to the 
c™P, which 



where A'^^p, B*^^p correspond to the exphcit part and A 



imp 



imphcit part. The methods share the same abscissa c 



cxp 



makes (2.11) internally consistent |29i. Definition 2]. The methods also share 
the same coefficient matrices U^^p = U^'^p = U and V^'^p = V'°'p = V. 
The coefficients q^^'' 



q™'' in (2.5) can be different, which means that the 



implicit and explicit components use different initialization and termination 



procedures. An IMEX-GLM (2.11) is a special case of a partitioned GLM [291 



Definition 1]; while in (2.11) the right hand side is split in two components, 
stiff and nonstiff, a partitioned GLM allows for splitting in an arbitrary 
number of components. 

It has been shown in [291 Theorem 2] that an internally consistent parti- 



tioned GLM (and, in particular, the IMEX GLM (2.11 )) has order p and stage 



order g G {p—l,p} if and only if each component method (A'^^p, B'^^p, U, V) 
and (A™p, B'™P, U, V) has order p and stage order q. We note that no 
additional "coupling" conditions are needed for the IMEX GLM (i.e., no or- 
der conditions that contain coefficients of both the implicit and the explicit 
schemes) . 



3 Extrapolation-based IMEX GLMs 

3.1 Method formulation 

In this section we derive the new extrapolation-based IMEX GLMs. Consider 
the following extrapolation formula depending on stage values Yj^ and Y^"" ' 



at two consecutive steps 



j-i 



fj = E".^/K"') + E/5^-^/W"^")' J- = i'2, . . . ,.. (3.1) 



n+l]^ 



k=l k=l 



Substituting fj in (3.1) for f{Yj ') in (2.1) leads to the proposed class 
of extrapolation-based IMEX GLMs. The simple example of IMEX method 
consisting of the explicit Euler method combined with the A-stable implicit 
^-method corresponding to ^ > 1/2 is presented in |18j . 



Substituting (3.1) into (2.1) leads to 



i j-1 



j=i fc=i j=i fc=i 

+ hj2^,,g{Y}-^'^)+j2u^df^ . = 1,2,...,., 



s s 



i=i 



,[n+l] 



i-1 



/^EE^^.«^-^/(^i'^') + ^EE^^./5^-^'/(^i"^") 



j=i fc=i j=i fc=i 

+ /^EMK-'^-^'VE^?/]"^' ^ = l,2,...,r, 

n = 0, 1,...,A^ — 1. Changing the order of summation in the double sums 
above and then interchanging the indices j and k we obtain 



i— 1 i 



[n+1] 



i=i k=i j=i k=j+i 

i r 

+ ^ E "-^^^ (^j"^^') + E ^^J'^i"^ ' ^ = 1> 2, . . . , s. 






i=i 



s— 1 s 



[n+1] 



= /^EE^^^^^^/ft-^'^V^E E &^^/3^./(^r'^) 



These relations lead to IMEX GLMs of the form 

+ /.^a,,^(r;"+^])+E«4"'' ^ = 1,2,...,., 

't i-l (3-2) 

i=i 3=1 

i=i i=i 

n = 0,l,...,A^ — 1, where the coefficients ajj, a*^-, 6jj, and 6*j- are defined by 

i i s s 

dij = 2_^'^ikOikj, O*,' = /_^ O-ikf^kj, bij = 2_^^ikOikj, Kj = 2_^ hkf^kj- 
k=l k=j+l fc=l k=j+l 

In matrix notation 

A = [a^J] G M^^^ A* = [ay G M^^^ B = [%] G M''^^ B* = [b*^] G M''^^ 

with 

A = Aa, A* = A/3, B = Ba, B* = B/3, 

where a = [aij] G M''^'', (3 = [f3ij] G M'*^''. Observe that the matrix A* is 
strictly lower triangular and that the last column of the matrix B* is zero. 
In matrix notation the extrapolation-based IMEX-GLM is defined by: 

f["+^] = /i(A®i)/(r["i) +/i(A*®i)/(r["+^]) 

+h{A ® I)(7(F["+^1) + (U I)i/["1, (3.3) 



y 



n+l] 



h{B ® I)/(Ff"]) + h{B* ® I)/(f["+']) 



+/i(B ® I)/ (r ["+!]) + (v®i)i/M, 

n = 0,l,...,Ar-l. 



The explicit part of (3.3), obtained for g{y) = 0, can be represented as a 
single GLM extended over two steps from t„_i to t„ and t„ to t„+i, as follows 



" yH " 




' 





I " 




y[n+l] 
y[n+l] 


= 


A 


A* 


U 




A 


A* 


U 


y[n+l] 




B 


B* 


V 





/(yN) 



y[nj 



(3.4) 



eV,c^V. 



The abscissa vector is c'^^p = [(c 

Similarly, the implicit part of the IMEX scheme (3.2) corresponding to 
f{y) = assumes the form 



yH 




" 


I " 




y[n+l] 
y[n+l] 


= 


A 


U 




A 


U 


^[n+1] 




B 


V 





^(yN) 

yN 

yH 



(3.5) 



This method has the order and stage order of the underlying GLM (2.1) 
since it is the same method. The abscissa vector is c'™p = [(c — e)^, c^^ 



and therefore the method (3.3) is internally consistent 



3.2 Construction of the interpolant 

We define the local discretization errors 77 (t„ + Cjh) of the extrapolation 



formula (3.1) by the relation 



f{yitn + Cjh)) = ^ajkf{yitn~i + Ckh)) + 



i-i 



fc=i 



(3.6) 



+ ^/^jkf{y{tn + Ckh)) + ri{tn + Cjh), 



fc=l 



j = 1,2, ... ,s. Letting (^(t) = f{y{t)) the relation (3.6) can be written in 
the form 

s j-l 

r]{tn + Cjh) = ip{tn + Cjh) - ^ ajk^[tn + (Cfc - l)h) - ^ (3jk^{tn + Ckh), 



k=l 



fc=l 



j = 1, 2, . . . , s. Expanding ip(tn + Cjh), ip{tn + (c^ - l)h), and ipit^ + Ckh) 
into Taylor series around t„ we obtain 



v{t^ + c,h) = x: (I - E «..^r^ - E /3.4 ) ^'^^'^^'") + ^^^'^')- 

1=0 ^ ■ fc=l 



fc=l 



Assuming that the extrapolation procedure given by (3.1) has order p, i.e., 



ri{tn + Cjh) = 0{hP), leads to the following system of equations for the inter- 
polation coefficients: 



i-i 



Ea,fc(cfe-lf = 4-E/3,,ci, ^ = 0,l,...,p-l, j = l,2,...,s. (3.7) 
fc=i fc=i 



In matrix notation we have 



0,l,...,p-l. 



(3.8) 



3.3 Stage and order conditions 



Lemma 3.1 Assume that the underlying GLM (2.4) has order p and stage 



order q = p or q = p — 1, and that the interpolation formula (3.1) has order 



p (3.8). Then the explicit method (3.4) has order p and stage order q 



Proof. 



The method 


(3.4) 


has the coefficients 








A^^P = 





, B^^P = 


A A*" 


, c'^'P = 


c — e 




A A* 


. 


B B* 




c 


where e = [1, . . . , 1]^ e W, 




U'=xp = 


" I " 

u 


, V'^p = 


" U " 
V 


; 


and the vectors 




qr = 


e 


cxp 

; q* = 


r(c-e)n 


, i = l,...,p. 






[ qo 






[ q* J 









We verify directly that the extrapolation-based explicit method (3.4) has 



stage order q, i.e., it satisfies equations (2.9) 





c^^ - A; A (« (c - e)'=-^ + (3 c^-^) - k\ Uqfc 






{from (3.8)} 







A; = 0,l,...,g {from (2.9)}. 



From (2.10) we verify that the method (3.4) has order p 



k\ 



Yl T\ ^"-' - ^ B"'' (c<=^P)^-^ - k\ V^^P q^,^P 



^=0 



£=0 



k\ 







A;! 



Uqfc 
Vqfc 



E^ ^'-"'ir -feAc^-^-fc!Uq, 



i=Q 



i\ {k-l)\ 

Y.^^k-e-kBc^-'-kWqk 



e=o 



{from (3.8)} 



k = 0,l,...,p {from (2.9), (2.10), and (3.9)} . 






For the first component we have used the fact that 



c^=((c-e) + e) =Em73w(^-^) 



k-e 



e=o 



k\{i-k)\ 



(3.9) 



D D 



To analyze the order and stage order of IMEX GLMs (3.2 ) we will impose 



some conditions on the local discretization errors of the internal and external 



10 



stages of the underlying GLM (2.1 ) and on the accuracy of the extrapolation 



procedure (3.1) 



We have the following result. 



Theorem 3.1 Assume that the underlying GLM (2.1) has order p and stage 
order q=porq = p — l, and that the extrapolation procedure (3.1) has order 
p. Then the IMEX GLM (3.2) has order p and stage order q = p or q = p — 1. 



Proof. 

The method gs) is an IMEX GLM of the form ( [2ll| ) pHJ Definition 
1]. The explicit (3.4) and implicit (3.5) components have the same order p 



and stage order q, and they share the same abscissa vector and the same 
coefficients 



,exp 



,imp 



u 



imp 



■yexp Y"'™P = \^^P 



The result follows directly from [281 Theorem 2]. 

D D 

3.4 Prothero-Robinson convergence of IMEX GLMs 



The extrapolation IMEX-GLM schemes (3.2) do not suffer from order re- 



duction phenomenon when applied to stiff systems of differential equations. 
Following |5l EH EH] we consider the Prothero-Robinson (PR) [21] test prob- 
lem of the form 




0(0), 



0(t))+0'(t), t>o, 



(3.10) 



where /i G C has a large and negative real part and 0(t) is a slowly varying 



function. The solution to (3.10) is y{t) = (j){t). The IMEX scheme (3.2) is 



said to be PR-convergent if the application of (3.2) to the equation (3.10) 



leads to the numerical solution y'"^ whose global error satisfies 



yt'^i-E^^^V'^Ht. 



fc=0 



0{h^) as /i — 7- and hjj, — )• — oo. 



We have the following result. 
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Theorem 3.2 Assume that the implicit GLM (2.1) has order p and stage 



order q = p — 1 or q = p, and that the extrapolation formula (3.1) has order 



p. Then the IMEX scheme (3.2) is PR-convergent with order niin(p, g) as 



/i — )■ 0, /i/i — )■ — oo, and hfi G Sj. Here, Sj is the stability region of the 



implicit GLM (3.5) 



Proof. The resuh foUows directly from 
of IMEX-GLMs. 

D D 



Theorem 3] on PR-convergence 



4 Stability analysis of IMEX GLMs 



To analyze stability properties of IMEX GLMs (3.2) we use the test equation 

y'it) = \oyit) + \Mt), t>0, (4.1) 



where Aq and Ai are complex parameters. Applying (3.3) to (4.1 ) and letting 
Zi = hXi, i = 0,1, we obtain 

(I - (^oA* + 2iA))r["+il = 2oAyW + UyM, 

-{zqB* + 2iB)F["+i] + ?/["+il = 2oBF["1 + VyM, 

n = 0,l,...,A^— 1. This is equivalent to the matrix recurrence relation 



y[n+l] 

yln+1] 



M{zo,zi] 



yH 



(4.2) 



where the stability matrix M(zo, ^i) is defined by 

M(^o,-2i) = 
with 



mii{zo,Zi) mi2{zo,Zi] 
m2i{zo,Zi) m22{zo,Zi] 



mu{zo, zi) = zo[l - {zqA* + ziA)) A, 
mi2izo,zi) = {I - {zqA* + ziA)y^V , 



12 



^21(^0, zi) = zo(b + {zqB* + ziB) (I - {zoA* + z^A)) ^A 
^22(^0, ^1) = V + {zoB* + z^B) (I - (20A* + z,A)y'v. 



We define also the stability function of the IMEX GLM (3.2) as a character- 
istic polynomial of the stability matrix M(zo, zi), i.e., 

p{w, zo, zi) = det {wl - M(zo, ^1)) • 

For Zq = the stability matrix M(0, Zi) and polynomial p(w, 0, zi) = w'^rj^w, zi) 
are those of the underlying GLM ( |2.1 ). For zi = we obtain M(zo,0), and 
it can be verified that this corresponds to the stability matrix of the explicit 



method (3.4). 



We say that the IMEX GLM ([3^ is stable for given ^o, -^i e C if all the 
roots Wi{zQ, zi), i = 1,2, . . . ,s + r, of the stability function p{w, zq, zi) are 
inside of the unit circle. As observed in [IHj in the context of IMEX 6'-methods 



of order one, a large region of absolute stability for the explicit method (3.4) 



and good stability properties (for example A- or L-stability) for the implicit 
method are not sufficient to guarantee desirable stability properties of the 



overall IMEX GLM (|3.2|). We have to investigate the stability properties of 

EHl. 



the combined as IMEX GLM 

In this paper we will be mainly interested in IMEX schemes which are 
A{a)- or A-stable with respect to the implicit part zi G C. To investigate 
such methods we consider, similarly as in [HI EH], the sets 



5„ 



zqEC : the IMEX GLM is stable for any 

2i G C : 7^(2l) < and |lm(2i)| < tan(a)|7e(zi) 



^a,y 



(4.3) 



(4.4) 



(4.5) 



Observe also that the region S^ is independent of a, and corresponds to the 



For fixed values of y G M we define also the sets 

Zo e C : the IMEX GLM is stable for fixed 

Zi = — |y|/tan(a) + iy 

It follows from the maximum principle that 



region of absolute stability of the explicit method (3.4). This region will be 
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denoted by Se- We have 

Sa C Se, (4.6) 

and we will look for IMEX GLMs for which the stability region Sa contains 



a large part of the stability region Se of the explicit method (3.4), for some 
a G (0, vr/2], preferably for a = 7i/2. 

The boundary dSa,y of the region Sa^y can be determined by the boundary 
locus method which computes the locus of the curve 

55„,j, = 1^0 G C : p{e'',Zo,-\y\/tiin{a)+ty)=0, ^ G [0, 2A;7r]}, 

where A; is a positive integer. 



2 
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W 1 ''^^^S^^^^^-^r---^^ 


^^ ^— — — N 
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-1.5 -1 

Re(z„) 



Figure 4.1: Points on the intersection of the ray yo = ttixq and dS.„/2,y 
for y = — 1.5, — 1.0, . . . , 1.5 (circles) and on the intersection of yo = itixq 
and dSTr/2 (square). This figure corresponds to IMEX GLM scheme with 
p = q = r = s = 2ioTm, = —1, A = 0.3, and /32i = 4.3 



We have also developed an algorithm to determine the boundary dSa of 
the stability region Sa- For fixed direction m corresponding to the ray 

yo = mxo, 

and for fixed Zi = — |y|/tan(a) + iy (or any Zi G C) we can compute the 
point Zq = Xo + iyo of intersection of the boundary dSa^y of Sa,y with the ray 
yo = mxo taking into account that such a point satisfies the condition 



max 

i=l,2,...,s+r 



Wi{zo,-\y\/tSiia{a) + iy) 



1. 
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This can be done by the bisection method which terminates when the con- 
dition 



max 

=l,2,...,s+r 



Wi[zQ,-\y\/iixn{a) + iy) 



1 



< toL 



(4.7) 



with accuracy tolerance tol is satisfied. We apply this method to the interval 



[xq, 0] with xq large enough so that the condition (4.7) is not satisfied for the 
first iteration of bisection method. This process leads to the definition of the 
function 

xo = f{m,a,y), 

where Xq corresponds to the points zq = Xo + iyo G dSa,y Then the boundary 
dSa of the region S^ can be determined by minimizing the negative value 
of this function for m G M and plotting the resulting points Zq = Xq + iyo- 



This algorithm is illustrated in Fig. |4.1[ where we have plotted points on 
the intersection of the ray yo = mxo and dS.„/2^y for y = —1.5, —1.0, . . . , 1.5 
(circles) and on the intersection oiyo = mxo and 95,^/2 (square). This figure 
corresponds to IMEX GLM scheme with p = q = r = s = 2 for m = —1, 
A = 0.3, and /32i = 4.3. This minimization can be accomplished using the 
subroutine fminsearch.m in Matlab applied to f{m, a, y) for fixed values of 
771 G M and a G [0, 7r/2] starting with appropriately chosen initial guesses for 
y. This process will be next applied to specific IMEX GLMs. 



Construction of IMEX GLMs with desir- 
able stability properties 



In this section we describe the construction of IMEX GLMs (3.2) up to the 
order p = 4 with large regions of absolute stability S with respect to the 
explicit part assuming that the implicit part is A- or L-stable. 

We will always start with implicit DIMSIM with p = q = r = s and the 
abscissa vector c given in advance. After computing the coefficient matri- 
ces A and V so that the resulting method has Runge-Kutta stability with 
the underlying Runge-Kutta formula which is A- or L-stable, the coefficient 
matrix B is computed from the relation 



B = Bo - ABi - VBa + VA. 



(5.1) 
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Here, Bq, Bi, and B2 are s x s matrices with the {i,j) elements given by 



l+Ci 



(t)j{x)dx 



:i + Q 



<H 



^j[Cj) 



(j)j{x)dx s 



i = 1,2, . . . ,s, compare Th. 5.1 in [6] or Th. 3.2.1 in [T9] . 



5.1 IMEX GLMs with p = q = r = s = l 

Consider the imphcit ^-method defined by 

" y[n+i] = /i^(/(F[«+i]) + ^(F["+il)) + yH 



Jn+ll 



/,(/(y["+i])+^(y[n+i]))+^N 



(5.2) 



n = 0, 1, . . . , A^ — 1. This method is A-stable for 6 G [1/2, 1] and L-stable for 
6 G (1/2, 1]. Consider also the extrapolation procedure 



f{Y 



[n+l]^ 



f{Y^' 



(5.3) 



n=l,2,...,A^-l. Substituting ^^ into ^^ we obtain IMEX 6'-method 
of the form 

(5.4) 

y[n+l] ^ /,(j(yN) ^^(y[n+l])) ^^N^ 

n = l,2,...,A^ — 1. We would like to point out that this variant of IMEX 
scheme is different from IMEX 6'- method considered in HHI. Observe that 



the method (5.4) requires a starting procedure to compute F^^l ~ j/(to + Oh) 
and y'^l ^ y{ti)- 



The method (5.2) can be represented by the abscissa c = 9, the parti- 
tioned matrix 



A 


U 


B 


V 



" 6 


1 


1 


1 



and qo = 1, qi = 0, and the explicit formula corresponding to g(y) = in 
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Q is GLM with c = [^ - 1, 9]' 









" 


1 " 


A 


U 


= 


e 


1 


B 


V 


e 


1 








. 1 


1 _ 



(5.5) 



andqo=[l,l]^, qi = [^-l,0]^ 



0.8 

0.6 
0.4 
0.2 









.^;<ei 


- 


^ 




^^ 


- 9=3/4 ^^^ 


'^X 




) 


^.^''^f 


>^e-2/3 


1 / e-1/2 


1 



-2.5 -2 -1.5 

Re(z„) 



Figure 5.1: Stability regions Se = Se{0) of explicit methods for 9 = 1/2, 
2/3, 3/4, 4/5, and 1 




Figure 5.2: Stability regions Sn/2,y{0), y = — 1.0, — 0.9, . . . , 1.0 (thin lines), 
<Sn/2{0) (shaded region), and Se{0) (thick line) for 9 = 3/4 



It can be verified that the stability matrix of IMEX scheme (5.4) takes 
the form 



M{zo,zi) 



l-dzi 



0zo 1 

zo 1 + (1 - 9)zi 
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Figure 5.3: Stability regions ST,/2,y{0), y = — 1.0, — 0.9, . . . , 1.0 (thin lines), 
<Stt/2{(^) (shaded region), and Se{0) (thick line) for 9 = 2/3 



To investigate stability properties of (5.4) it is more convenient to work with 



the polynomial obtained by multiplying the characteristic function p{w, zq, zi) 
of M(zo, zi) by a factor {l — 6zi)^. The resulting quadratic polynomial, which 
will be denoted by the same symbol p{w, Zq, Zi), takes the form 

p{w,zo,zi) = {i-ezifw^-{i-ezi){i+ezo+{i-e)zi)w-{i~e)zo{i-ezi). 



The stability polynomial of the explicit methods ( 5.5[ ) obtained by letting 
g{y) = in (5.4:) corresponds to 2:1 = and is given by 



p{w, Zq, 0) = W^ - {1 + 9zq)w - (1 - 6)Zo. 

The stability regions Se = Se{0) corresponding to this polynomial are plot- 



ted in Fig 5.1 for 6 = 1/2, 2/3, 3/4, 4/5 and 1. Observe that the stability 



region of the method (5.5) corresponding to = 1 is the unit disk 



Se = Se{1) = {zoeC: \zo + l\ <1}. 



We will investigate next the regions 5^/2 = 5^/2(6') defined by (4.3). We 



analyze first the case ^ = 1. It follows from Schur criterion applied to the 
polynomial p{w, Zq, iy) that zq G 5^/2(1) if and only if 

y^ - 2xo - Xq - ?/o > 

for any y eM.. This is equivalent to 

{xo + lf + yl<l or |zo + l|<l 



and it follows that in this case St,/2{1) = Se{1)- For 9 G (1/2, 1) St,/2{0) 
is no longer equal to Se{0), but St,/2{0) contains some part of Se{0). This 
is illustrated in Fig. 5.2| for 9 = 3/4 and in Fig. 5.3 for 9 = 2/3. In these 



figures we have plotted the regions S.„/2,y = 5^/2,j/(6') defined by (4.4) for 
y = — 1.0, — 0.9, . . . , 1.0 (thin lines), the regions 5^/2 = i5^/2(^) defined by 
(4.3) (shaded regions), and the regions Se = <Se{9) corresponding to explicit 
methods (5.5) for 9 = 3/4 and 9 = 2/3. We can see that for these values of 



9 the sets St,/2{9) contain quite large parts of Se{9). These figures illustrate 
also the relation (4.5) and the inclusion (4.6). For 9 = 1/2 it follows from 
Schur criterion that zq G iS7r/2(l/2) if and only if 



X, 



+ yl<4: and x^y'^ - Ayyo - Xoixo + 2^ - y^A + xq) > 



for any y G M. This implies that 

xo > and (2 + xo)^(xo + ^q) < 
and it follows that 5^/2(1/2) is empty. 

5.2 IMEX GLMs with p = q = r = s = 2 

Consider the implicit DIMSIM with c = [0, 1]"'", the coefficient matrices given 

by 



A 
B 



U 
V 



A 

2 



1+2A 
8A^+12A^-2A+5 

4(2A+1) 
8A3+20A2-2A+3 




A 


1 
1 


1-4A2 
4 


l + A 1- 



-8A3-12A2+10A-1 



4(2A+1) 

and the vectors qo, qi, and q2 equal to 



4(2A+1) 



o + A 7s — A 



qo 



qi 



-A 

-2A^+A-l 
2A+1 



q2 





1-2A 



It was demonstrated in [19] that this method has order p = 2 and stage 
order q = 2. Moreover, this method is A-stable if A > 1/4 and L-stable for 

A = (2 ± V2)/2. 
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The coefficients ajk of the extrapolation formula (3.1) of order p = 2 



computed from the system (3.7) corresponding to p = s = 2 take the form 



a 



ail ai2 

^21 ^22 



1 

-1 2-/321 



and the matrices A, A*, B, and B* appearing in the representation of IMEX 
GLM (3.2), and the corresponding explicit scheme (3.3) or (3.4), are gi^ 
by 



nven 



A 



B 



A 

, 2+(2-fei)A+2(2-fei)A^ 
^ 1+2A 





/32iA 



4A2-1 7-fei+2(l-fei)A+4(l+/32i)A2-8(l-fei)A3 

4 4(1+2A) 

1-10A+12A2+8A3 l+/92i+2(9-5/32i)A-4(l-3/32i)A^-8(l-fei)A=' 

4(1+2A) 4(1+2A) 



B* 



/32i(1-4A2) ^ 

4 ^ 

-fei(l-10A+12A^+8A3) 

4(1+2A) 













- 


area of S^ 

arpa nf R 


" 




1 ' 



Figure 5.4: Areas of stability regions Se 
a = 71 /2 and (32i e [0, 8] 



Se{(32i) and Sa = Sa{f32i) for 



To investigate the stability properties of the resulting IMEX scheme we 
will work with the stability polynomial p{w,Zo,zi) obtained by multiplying 
stability function of the method by a factor (1 — Xzi)"^. It can be verified 
that this polynomial takes the form 

p{w,zo,zi) = w{{l - Xzifw^ -p3{zo,zi)w'^ + p2{zo, zi)w - pi{zo, Zi)) , 

20 




Figure 5.5: Stability regions S.„/2,y, y = —2.0, —1.8, . . . , 2.0 (thin lines), 5-^/2 
(shaded region), and Se (thick line) for A = (2 — v^)/2 and /32i ~ 4.64 




Figure 5.6: Contour plots of the area of stability region St^/2 of IMEX GLMs 
for p = q = r = s = 2 



with the coefficients pi{zo, zi), P2{zq, zi), and Ps^zq, zi) which depend also on 
A and /32i- These coefficients are given by 



Pi{zo,zi) 



4-/32i-2(4-/32i)A + 4/32iA2-8/32iA3 1 - 4A + 2A 



4(1 + 2A) 
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Zo + 



^0> 




Figure 5.7: Stability regions S.„/2,y, y = —2.0, —1.8, . . . , 2.0 (thin lines), iS,r/2 
(shaded region), and Se (thick line) for A ~ 0.29 and /32i ~ 4.64 



P2{Zq,Zi) 



2(l-A + (2 



/32i)A2 - 2/32iA3) /321 
Zq H 



4/32iA + 2(l + /32i)A2 



2-/3: 
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1 + 2A 

4(2 - /32i)A + 2(2 



/32i)A^ 



-^O^^l) 



P3l2;o,^iJ 



1 + 



8 + /321 + 2(4 - /32i)A + 4(4 - 3/32i)A2 - 8/32iA=^ 



+ 



P2i\^zl - (2 



4(1 + 2A) 

P2i)zqz^ + (1 - 2A)2i 



-2;o 



1 - 4A + 2A2 



The underlying implicit GLM is A- and L-stable for A = (2 ± a/2)/2, 
compare [19], and we choose A = (2 — v^)/2 since this value leads to explicit 
methods and IMEX schemes with larger regions of stability Se and S.,^j2 than 
those corresponding to A = (2 + v2)/2. We have plotted in Fig 5.4 the area 
of the stability region Se = Se{P2i) of the explicit method (corresponding 
to zi = 0) and the area of the stability region St,/2 = 5^/2 (/32i) of the IMEX 
scheme for /32i G [0,8]. It can be verified that the explicit formula attains 
the maximal area of Se, approximately equal to 7.15 for f32i ~ 4.56, and the 
IMEX scheme attains the maximal area of iS7r/2, approximately equal to 5.75 
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4.64. 



for/3; 

On Fig 5.5 we have plotted stability regions ST,/2,y for y = —2.0, —1.8, . . . , 2.0 
(thin lines), stability region 5^/2 (shaded region), and stability region Se 
(thick line). We can see that iS^/2 contains a significant part of Se- 



We have also displayed on Fig. |5.6| contour plots of the area of stability 
region S^/2 of IMEX methods for /32i G [3,6] and A e [0.25,0.35]. This 
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area attains its maximum value approximately equal to 5.83 for /32i ~ 4.59 



and A ~ 0.29. This point is marked by the symbol 'x' on Fig. 5.6 On 



Fig 5.7 we have plotted stability regions 15,^/2,2; for y = —2.0, —1.8, . . . ,2.0 
(thin lines), stability region St,/2 (shaded region), and stability region Se 
(thick line) corresponding to these values of /32i and A. We can see again 
that i57r/2 contains a significant part of Se- 



5.3 IMEX GLMs with p = q = r = s = 3 

Let A ~ 0.43586652 be a root of the cubic polynomial 



y,(A) = A^ - 3A2 + -A - 



6' 



and consider the implicit DIMSIM with c = [0, 1/2, 1]"^, the coefficient matrix 
A given by 

0.43586652 

0.25051488 0.43586652 

-1.2115943 1.0012746 0.43586652 

the rank one coefficient matrix V = ev^, where e = [1, 1, 1]-^ and 

0.55209096 0.73485666 -0.28694762 
and the vectors qo, qi, q2, and qs equal to qo = e. 



qi 



-0.43586652 -0.18638140 0.77445315 



q2 
qs 



iT 



-0.092933261 -0.43650382 



-0.033649982 -0.17642592 



Computing the coefficient matrix B from the relation (5.1) leads to the 



method of order p = 3 and stage order q = 3. It was demonstrated in 
[19] that the resulting method is A- and L-stable. 
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The coefficients ajk of the extrapolation formula (3.1) of order p = 3 



computed from the system (3.7) corresponding to p = s = 3 take the form 



a 



"11 "12 "13 
a2l "22 «23 
asi "32 «33 




1 

3-/332 3/3: 



32 



1 

3 3-/321 

-8 6-/331-3/332 



To investigate the stability properties of the resulting IMEX scheme we 
will work with the stability polynomial p{w, zo,zi) obtained by multiplying 
stability function of the method by a factor {1 — XziY, where A is the diagonal 
element of the matrix A. It can be verified that this polynomial takes the 
form 

p{w,Zo,Zi) = {I- Xzifw^ -p;^{zo,Zi)w^ +p4{zo,Zi)w'^ -p3{zo,Zi)w^ 
+ P2izo,zi)w'^ - pi{zo, Zi)w + Po{zo, Zi), 

with the coefficients ^0(2:0, ^1), Pi{zo, zi), p2{zo, Zi), P3{zo, zi), Pa{zo, zi), and 
^5(2:0, Zi) which are polynomials of degree less than or equal to 3 with respect 
to Zq and Zi. These coefficients depend also on /32i, (3^1, and /332. 




Figure 5.8: Stability regions S.„i2,y, y = —2.0, —1.8, . . . , 2.0 (thin lines), 57r/2 
(shaded region), and Se (thick line) for /32i ~ 1.13, P-n ~ 1.45, and /332 ~ 
-0.158 



We have performed a computer search in the parameter space /32i, /33i, 
and P32 looking first for methods for which the stability region Se of the 
explicit method is maximal. This corresponds to the parameter values /32i ~ 
1.13, /331 K. 1.45, /332 ~ —0.158, for which the area of Se is approximately 
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Figure 5.9: Stability regions S.„/4^^y, y = —2.0, —1.8, . . . , 2.0 (thin lines), 5^/4 
(shaded region), and Se (thick line) for /32i ~ 1.13, f3si ~ 1.45, and f3s2 ~ 
-0.158 




Figure 5.10: Stability regions S.„i2,y, y = —2.0, —1.8, ... ,2.0 (thin lines), 
5,r/2 (shaded region), and Se (thick line) for (32i ~ 1.39, (3-^1 ~ —0.146, and 
/332 ~ 1.24 



equal to 3.54. The stability region Se of the resulting method is plotted 
on Fig. 5.8 by a thick line. We have also plotted stability regions 15,^/2,2/ 
for y = —2.0, —1.8, . . . ,2.0 (thin lines) and the stability region 5^/2 (shaded 
region) of the corresponding IMEX scheme. We can see that this region 5^/2 
is substantially smaller than the region Se, the area of i5,r/2 is approximately 
equal to 0.39. However, we can obtain larger regions Sa for values of a 
smaller than 7r/2, i.e., if we relax the requirement that the implicit part of 
IMEX scheme is A-stable and require instead yl(a)-stability for a < 7r/2. 
This is illustrated on Fig. 5.9, where we have plotted again the stability 
region Se of the explicit method (thick line), stability regions of Sa,y for 
y = —2.0, —1.8, . . . , 2.0 (thin lines), and stability region Sa (shaded region) 
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Figure 5.11: Stability regions STr/4^y, y = —2.0, —1.8, ... ,2.0 (thin lines), 
5^/4 (shaded region), and Se (thick line) for /32i ~ 1.25, /^ai ~ 1.62, and 



A 



>7r/4 
32 



0.0555 



for a = 7r/4. The area of this region is approximately equal to 1.91. 

We have also performed a computer search looking for methods for which 
stability regions Sa are maximal for some fixed values of a. For a = 7r/2 
this corresponds to the parameter values /32i ~ 1.39, (3^1 ~ —0.146, and 
P32 ~ 1-24 for which the area of (5,^/2 is approximately equal to 0.50. For 
a = 7r/4 this corresponds to the parameter values /32i ~ 1.25, P^i ~ 1.62, and 
(3^2 ~ 0.00555 for which the area of iS,r/4 is approximately equal to 2.80. We 
have plotted on Fig. |5.10| and Fig. 5.11 the stability regions Se of the resulting 
explicit methods (thick lines), the regions Sa,y for y = — 2.0, — 1.8, . . . , 2.0 
(thin lines) and stability regions Sa of IMEX schemes (shaded regions) for 
a = tt/2 and a = tc/A. 



5.4 IMEX GLMs with p = q = r = s = 4 

Let A ~ 0.57281606 be a root of the polynomial 

^(A) = A^-4A^ + 3A^-^A + ^, 
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and consider the implicit DIMSIM with c = [0, 1/3, 2/3, 1]^, the coefficient 
matrix A given by 



A 



0.57281606 

0.15022075 

0.59515808 

1.7717286 





0.57281606 

-0.26632807 0.57281606 



-1.64234444 0.39147320 0.57281606 
the rank one coefficient matrix V = ev"^, where e = [1, 1, 1, 1]-^ and 

15.615037 -46.967269 41.290082 -8.9378502 
and the vectors qo, qi, q2, qs, and q4 equal to qo = e. 



iT 



qi 



q2 

qs = 

q4 = 



-0.57281606 -0.38970348 -0.23497940 

-0.13538313 -0.070879128 0.21364995 



0.093673420 



-0.025650275 -0.063113738 -0.11549405 



-0.0030214983 -0.018412760 -0.062996758 



Computing the coefficient matrix B from the relation (5.1) leads to the 



method of order p = 4 and stage order g = 4. It was demonstrated in 
that the resulting method is A- and L-stable. 



The coefficients ajk of the extrapolation formula (3.1) of order p = 4 



computed from the system (|3.7|) corresponding to p = s = 4 take the form 

1 



a 



with 





-1 4 -6 4-/321 

1332-4: 15-4/332 2(3/332-10) 10-/331-4/332 

a^i 0^42 0^43 044 



^41 = /342 + 4/343 - 10, 042 = 36 - 4/342 - 15/343, 

043 = 6/342 + 20/343 - 45, a44 = 20 - /34i - 4/342 - IO/343. 
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Figure 5.12: Stability regions 5,^/4,2;, y = —2.0, —1.8, ... ,2.0 (thin lines), 
57r/4 (shaded region), and Se (thick line) for /32i ~ 0.0645, /Ssi ~ —0.351, 
and 1332 ^ 0.272, /34i ~ -2.82, 13^2 ~ 3.47, /343 ~ -1-05 




Figure 5.13: Stability regions Sa for a = 0, 7r/12, 7r/6, 7r/4, 7r/3, 57r/12, 
and 7r/2 for /^si ~ 0.0645, /^si ^ -0.351, and /332 ^ 0.272, /34i ~ -2.82, 
/342 ~ 3.47, /343 ~ -1.05 



To investigate the stability properties of the resulting IMEX scheme we 
will work with the stability polynomial p{w,Zo,zi) obtained by multiplying 
stability function of the method by a factor {1 — XziY, where A is the diagonal 
element of the matrix A. It can be verified that this polynomial takes the 
form 

p{w,zo,zi) = {1- Xzi^w^ -p7{zo,zi)w'^ +pq{zo,zi)w^ -p5{zo,zi)w^ 

+ P4.{zo, zi)w^ - psizo, zi)w^ + p2{zo, Zi)w'^ - pi{zo, Zi)w + po{zo, Zi), 

with coefficients po{zo, Zi), pi{zo, Zi), ^2(^0, ^1), ^3(^0, ^1), P4(^o, ^1), ^5(^0, ^1), 
Pq{zo,Zi), and pj{zq,Zi) which are polynomials of degree less than or equal 
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Figure 5.14: Stability regions S.„i2,y, y = —2.0, —1.8, ... ,2.0 (thin lines), 
5^/2 (shaded region), and Se (thick hne) for /32i ~ —0.00516, /S^i ^ —0.939, 
/332 ~ 1.18, /341 ^ -1.71, /342 ~ 2.07, and (3^3 ~ 0.32 




Figure 5.15: Stability regions S.„/4,^y, y = —2.0, —1.8, ... ,2.0 (thin lines), 
5^/4 (shaded region), and Se (thick line) for /32i ~ 0.0964, P^i fti —0.278, 
(332 ~ 0.464, /341 ^ -1.63, /342 ~ 2.73, and (3^3 ^ -0.678 



to 4 with respect to 2^0 and zi. These coefficients depend also on (32i, /Ssi, 
(332, (341, (3 42, and (343- 

We have performed a computer search in the parameter space (32i, /Ssi, 
P32, (341, (342, and (3^3 looking first for methods for which the stability region 
Se of the explicit method is maximal. This corresponds to the parameter 
values (321 ^ 0.0625, /Jgi ^ -0.355, (332 ^ 0.272, (3^i ^ -2.84, (342 ^ 3.49, 
(3^3 ^ —1.06, for which the area of Se is approximately equal to 2.82. The 



stability region Se of the resulting method is plotted on Fig. |5.12| by a thick 
line. We have also plotted stability regions 5^/4,^ for y = —2.0, —1.8, . . . , 2.0 
(thin lines) and the stability region 5^/4 (shaded region) of the corresponding 
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IMEX scheme. The area of iS,r/4 is approximately equal to 0.32. We have 



also plotted on Fig. 5.13 stability regions Sa for a = 0, vr/12, n/Q, 7r/4, 
7r/3, 57r/12, and 7r/2 corresponding to the same values of (3ij. We can see in 
particular that the region 5,r/2 is quite small, its area is approximately equal 
to 0.0069. 



As in Section 5.3 we have also performed a computer search looking 
directly for methods for which stability regions Sa are maximal for some 
fixed values of a. For a = 7r/2 this corresponds to the parameter values 
1321 ~ -0.00516, /331 ^ -0.939, /332 ^ 1-18, (3^1 ^ -1.71, /342 ^ 2.07, and 
/343 ~ 0.32, for which the area of iS7r/2 is approximately equal to 0.16. For a = 
7r/4 this corresponds to the parameter values /32i ~ 0.0964, (3^1 ^ —0.278, 
(332 ~ 0.464, /341 ^ -1.63, ^42 ~ 2.73, and ^43 ~ -0.678, for which the 



area of 5,^/4 is approximately equal to 0.65. We have plotted on Fig. 5.14 



and Fig. 5.15 the stability regions Se of the resulting explicit methods (thick 
lines), the regions Sa,y for y = —2.0, —1.8, . . . , 2.0 (thin lines) and stability 
regions S^ of IMEX schemes (shaded regions) for a = tt/2 and a = 7r/4. 



6 Numerical experiments 

The extrapolation-based IMEX GLMs constructed in Section |5] have been 
implemented in Matlab. The required starting values y^^^ and Y'^ were com- 
puted by finite difference approximations from solutions obtained with the 
Matlab routine odel5s. 

The test problem is the two dimensional shallow-water equations system 
|21j . which approximates a thin layer of fluid inside a shallow basin: 

I" + !<""' + !<'"') = " 

|(.ft) + |;(„,^ft + ljft=)+|(»A) = (6,1) 

Here h(t, x, y) is the fluid layer thickness, u{t, x, y) and f (t, x, y) are the com- 
ponents of the velocity field, and g denotes the gravitational acceleration. 
The spatial domain is fi = [—3, 3]^ (spatial units), and the integration win- 
dow is to = < t < tf = 10 (time units). We use reflective boundary 



30 



conditions and the initial conditions at to = 

u{to, X, y) = v{to, x,y)=0, h{to, x,y) = l + e-ll(-.^)-(-i''=2)lli , (6.2) 

with the Gaussian height profile ci = 1/3 and C2 = 2/3. 

A second order Lax-Wendroff finite difference scheme is used for space 
discretization, resulting in a semi-discrete ODE system of the form 

^Uit) = F{U)= Fu{U) ■ f/(t) + (F(f/) - Fu{U) ■ ^(t)) , (6.3) 

9iU) f(U) 

where U{t) is a combined column vector of the discretized state variables 
(/i, m/i, vh), and Fu = dF/dU be the Jacobian of right hand side. We consider 
a sphtting of equation ( 6.3[ ) into the linear part g{U), considered stiff, and 



the nonlinear part f{U), considered non-stiff. The linear stiff part is treated 
implicitly, and the non-stiff part is treated explicitly. 

We compare the numerical results for the solution at the final time with a 
reference solution computed by the Matlab function ode 15s with very tight 
tolerances atol = rtol = 10"^"^. The errors are measured in £2 norms. The 



error diagram against the time step size is presented in Fig. 6.1 The observed 



orders for all the methods tested match the theoretical predictions. 

7 Concluding remarks 

General linear methods offer an excellent framework for the construction of 
implicit-explicit schemes. In this paper we develop a new extrapolation-based 
approach for the construction of practical IMEX GLM schemes of high order 
and high stage order. The accuracy, linear stability, and Prothero-Robinson 
convergence are analyzed. These schemes are particularly attractive for solv- 
ing stiff problems, where other multistage methods may suffer from order 
reduction. 

The extrapolation-based mechanism offers a systematic approach for con- 
structing IMEX GLM schemes. The construction starts with the selection of 
an implicit component method with suitable stability and order properties. 
The explicit component is then obtained though an optimization procedure 
that maximized the combined stability region of the pair. We apply this 
methodology to construct IMEX pairs of orders one to four. 
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Figure 6.1: Error vs. time step size for several extrapolation-based IMEX- 
GLMs applied to the shallow water equations test problem. 



Future work is planned to extend the extrapolation idea to construct 
other types of partitioned GLMs, including parallel time integrators, and 
asynchronous pairs of methods that do not share the same abscissae. 
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